### load packages and data
library(pacman)
p_load(tidyverse, tidycensus, sf)

rootdir <- getwd()
setwd(paste(rootdir, "/Data/ACS/", sep=""))
## connect to census API (you need a personal API key)
source("../../census_API.R")

# variable names
v8<-load_variables(2018, "acs5", cache = TRUE)

NYC_counties <- c("36005", "36047", "36061", "36081", "36085")

# pull acs data by race for years 2014 - 2018 (5 year survey)
## H variables are white non hispanic
## B variables are black but include black hispanic
## D variables are asian
## I variables are all hispanic

for (i in c(2014:2018)) {
        acs_white <- get_acs(state = "NY", geography = "tract", year = i,
                         variables = c(pop_white_nonhisp = "B03002_003",
                                       male_und_5 = "B01001H_003",
                                       male_9 = "B01001H_004",
                                       male_14 = "B01001H_005",
                                       male_17 = "B01001H_006",
                                       male_19 = "B01001H_007",
                                       male_24 = "B01001H_008",
                                       male_29 = "B01001H_009",
                                       male_34 = "B01001H_010",
                                       male_44 = "B01001H_011",
                                       male_54 = "B01001H_012",
                                       male_64 = "B01001H_013",
                                       male_74 = "B01001H_014",
                                       male_84 = "B01001H_015",
                                       male_85up = "B01001H_016",
                                       fem_und_5 = "B01001H_018",
                                       fem_9 = "B01001H_019",
                                       fem_14 = "B01001H_020",
                                       fem_17 = "B01001H_021",
                                       fem_19 = "B01001H_022",
                                       fem_24 = "B01001H_023",
                                       fem_29 = "B01001H_024",
                                       fem_34 = "B01001H_025",
                                       fem_44 = "B01001H_026",
                                       fem_54 = "B01001H_027",
                                       fem_64 = "B01001H_028",
                                       fem_74 = "B01001H_029",
                                       fem_84 = "B01001H_030",
                                       fem_85up = "B01001H_031",
                                       no_ins_18und = "C27001H_004",
                                       no_ins_19_64 = "C27001H_007",
                                       no_ins_65_up = "C27001H_010",
                                       households = "B11001H_001",
                                       family_hhs = "B11001H_002",
                                       married_couple_hhs = "B11001H_003",
                                       living_alone = "B11001H_008", 
                                       multigenerational = "B10051H_001",
                                       inc_und_10k = "B19001H_002",
                                       inc_15k = "B19001H_003",
                                       inc_20k = "B19001H_004",
                                       inc_25k = "B19001H_005",
                                       inc_30k = "B19001H_006",
                                       inc_35k = "B19001H_007",
                                       inc_40k = "B19001H_008",
                                       inc_45k = "B19001H_009",
                                       inc_50k = "B19001H_010",
                                       inc_60k = "B19001H_011",
                                       inc_75k = "B19001H_012",
                                       inc_100k = "B19001H_013",
                                       inc_125k = "B19001H_014",
                                       inc_150k = "B19001H_015",
                                       inc_200k = "B19001H_016",
                                       inc_over_200k = "B19001H_017",
                                       agg_income = "B19025H_001",
                                       owner_occupied = "B25003H_002",
                                       pplxroom_morethan1 = "B25014H_003",
                                       foodstamps = "B22005H_002",
                                       tot_educ_25up = "C15002H_001",
                                       lesshigh_m = "C15002H_003",
                                       lesshigh_f = "C15002H_008",
                                       highschool_m = "C15002H_004",
                                       highschool_f = "C15002H_009",
                                       somecollege_m = "C15002H_005",
                                       somecollege_f = "C15002H_010",
                                       bahigher_m = "C15002H_006",
                                       bahigher_f = "C15002H_011",
                                       occupation_tot_16_up = "C24010H_001",
                                       service_occ_male = "C24010H_004",
                                       production_occ_male = "C24010H_007",
                                       service_occ_female = "C24010H_010",
                                       production_occ_female = "C24010H_013",
                                       tot_pop_occupied_housing = "B25032H_001",
                                       units_oned = "B25032H_002",
                                       units_onea = "B25032H_003",
                                       units_two = "B25032H_004",
                                       units_three_four = "B25032H_005"))
        acs_black <- get_acs(state = "NY", geography = "tract", year = i,
                               variables = c(pop_black = "B02001_003", 
                                             male_und_5 = "B01001B_003",
                                             male_9 = "B01001B_004",
                                             male_14 = "B01001B_005",
                                             male_17 = "B01001B_006",
                                             male_19 = "B01001B_007",
                                             male_24 = "B01001B_008",
                                             male_29 = "B01001B_009",
                                             male_34 = "B01001B_010",
                                             male_44 = "B01001B_011",
                                             male_54 = "B01001B_012",
                                             male_64 = "B01001B_013",
                                             male_74 = "B01001B_014",
                                             male_84 = "B01001B_015",
                                             male_85up = "B01001B_016",
                                             fem_und_5 = "B01001B_018",
                                             fem_9 = "B01001B_019",
                                             fem_14 = "B01001B_020",
                                             fem_17 = "B01001B_021",
                                             fem_19 = "B01001B_022",
                                             fem_24 = "B01001B_023",
                                             fem_29 = "B01001B_024",
                                             fem_34 = "B01001B_025",
                                             fem_44 = "B01001B_026",
                                             fem_54 = "B01001B_027",
                                             fem_64 = "B01001B_028",
                                             fem_74 = "B01001B_029",
                                             fem_84 = "B01001B_030",
                                             fem_85up = "B01001B_031",
                                             no_ins_18und = "C27001B_004",
                                             no_ins_19_64 = "C27001B_007",
                                             no_ins_65_up = "C27001B_010",
                                             households = "B11001B_001",
                                             family_hhs = "B11001B_002",
                                             married_couple_hhs = "B11001B_003",
                                             living_alone = "B11001B_008", 
                                             multigenerational = "B10051B_001",
                                             inc_und_10k = "B19001B_002",
                                             inc_15k = "B19001B_003",
                                             inc_20k = "B19001B_004",
                                             inc_25k = "B19001B_005",
                                             inc_30k = "B19001B_006",
                                             inc_35k = "B19001B_007",
                                             inc_40k = "B19001B_008",
                                             inc_45k = "B19001B_009",
                                             inc_50k = "B19001B_010",
                                             inc_60k = "B19001B_011",
                                             inc_75k = "B19001B_012",
                                             inc_100k = "B19001B_013",
                                             inc_125k = "B19001B_014",
                                             inc_150k = "B19001B_015",
                                             inc_200k = "B19001B_016",
                                             inc_over_200k = "B19001B_017",
                                             agg_income = "B19025B_001",
                                             owner_occupied = "B25003B_002",
                                             pplxroom_morethan1 = "B25014B_003",
                                             foodstamps = "B22005B_002",
                                             tot_educ_25up = "C15002B_001",
                                             lesshigh_m = "C15002B_003",
                                             lesshigh_f = "C15002B_008",
                                             highschool_m = "C15002B_004",
                                             highschool_f = "C15002B_009",
                                             somecollege_m = "C15002B_005",
                                             somecollege_f = "C15002B_010",
                                             bahigher_m = "C15002B_006",
                                             bahigher_f = "C15002B_011",
                                             occupation_tot_16_up = "C24010B_001",
                                             service_occ_male = "C24010B_004",
                                             production_occ_male = "C24010B_007",
                                             service_occ_female = "C24010B_010",
                                             production_occ_female = "C24010B_013",
                                             tot_pop_occupied_housing = "B25032B_001",
                                             units_oned = "B25032B_002",
                                             units_onea = "B25032B_003",
                                             units_two = "B25032B_004",
                                             units_three_four = "B25032B_005"))
        acs_asian <- get_acs(state = "NY", geography = "tract", year = i,
                             variables = c(pop_asian = "B02001_005", 
                                           male_und_5 = "B01001D_003",
                                           male_9 = "B01001D_004",
                                           male_14 = "B01001D_005",
                                           male_17 = "B01001D_006",
                                           male_19 = "B01001D_007",
                                           male_24 = "B01001D_008",
                                           male_29 = "B01001D_009",
                                           male_34 = "B01001D_010",
                                           male_44 = "B01001D_011",
                                           male_54 = "B01001D_012",
                                           male_64 = "B01001D_013",
                                           male_74 = "B01001D_014",
                                           male_84 = "B01001D_015",
                                           male_85up = "B01001D_016",
                                           fem_und_5 = "B01001D_018",
                                           fem_9 = "B01001D_019",
                                           fem_14 = "B01001D_020",
                                           fem_17 = "B01001D_021",
                                           fem_19 = "B01001D_022",
                                           fem_24 = "B01001D_023",
                                           fem_29 = "B01001D_024",
                                           fem_34 = "B01001D_025",
                                           fem_44 = "B01001D_026",
                                           fem_54 = "B01001D_027",
                                           fem_64 = "B01001D_028",
                                           fem_74 = "B01001D_029",
                                           fem_84 = "B01001D_030",
                                           fem_85up = "B01001D_031",
                                           no_ins_18und = "C27001D_004",
                                           no_ins_19_64 = "C27001D_007",
                                           no_ins_65_up = "C27001D_010",
                                           households = "B11001D_001",
                                           family_hhs = "B11001D_002",
                                           married_couple_hhs = "B11001D_003",
                                           living_alone = "B11001D_008", 
                                           multigenerational = "B10051D_001",
                                           inc_und_10k = "B19001D_002",
                                           inc_15k = "B19001D_003",
                                           inc_20k = "B19001D_004",
                                           inc_25k = "B19001D_005",
                                           inc_30k = "B19001D_006",
                                           inc_35k = "B19001D_007",
                                           inc_40k = "B19001D_008",
                                           inc_45k = "B19001D_009",
                                           inc_50k = "B19001D_010",
                                           inc_60k = "B19001D_011",
                                           inc_75k = "B19001D_012",
                                           inc_100k = "B19001D_013",
                                           inc_125k = "B19001D_014",
                                           inc_150k = "B19001D_015",
                                           inc_200k = "B19001D_016",
                                           agg_income = "B19025D_001",
                                           inc_over_200k = "B19001D_017",
                                           owner_occupied = "B25003D_002",
                                           pplxroom_morethan1 = "B25014D_003",
                                           foodstamps = "B22005D_002",
                                           tot_educ_25up = "C15002D_001",
                                           lesshigh_m = "C15002D_003",
                                           lesshigh_f = "C15002D_008",
                                           highschool_m = "C15002D_004",
                                           highschool_f = "C15002D_009",
                                           somecollege_m = "C15002D_005",
                                           somecollege_f = "C15002D_010",
                                           bahigher_m = "C15002D_006",
                                           bahigher_f = "C15002D_011",
                                           occupation_tot_16_up = "C24010D_001",
                                           service_occ_male = "C24010D_004",
                                           production_occ_male = "C24010D_007",
                                           service_occ_female = "C24010D_010",
                                           production_occ_female = "C24010D_013",
                                           tot_pop_occupied_housing = "B25032D_001",
                                           units_oned = "B25032D_002",
                                           units_onea = "B25032D_003",
                                           units_two = "B25032D_004",
                                           units_three_four = "B25032D_005"))
        acs_hisp <- get_acs(state = "NY", geography = "tract", year = i,
                            variables = c(pop_hisp = "B03002_012", 
                                          male_und_5 = "B01001I_003",
                                          male_9 = "B01001I_004",
                                          male_14 = "B01001I_005",
                                          male_17 = "B01001I_006",
                                          male_19 = "B01001I_007",
                                          male_24 = "B01001I_008",
                                          male_29 = "B01001I_009",
                                          male_34 = "B01001I_010",
                                          male_44 = "B01001I_011",
                                          male_54 = "B01001I_012",
                                          male_64 = "B01001I_013",
                                          male_74 = "B01001I_014",
                                          male_84 = "B01001I_015",
                                          male_85up = "B01001I_016",
                                          fem_und_5 = "B01001I_018",
                                          fem_9 = "B01001I_019",
                                          fem_14 = "B01001I_020",
                                          fem_17 = "B01001I_021",
                                          fem_19 = "B01001I_022",
                                          fem_24 = "B01001I_023",
                                          fem_29 = "B01001I_024",
                                          fem_34 = "B01001I_025",
                                          fem_44 = "B01001I_026",
                                          fem_54 = "B01001I_027",
                                          fem_64 = "B01001I_028",
                                          fem_74 = "B01001I_029",
                                          fem_84 = "B01001I_030",
                                          fem_85up = "B01001I_031",
                                          no_ins_18und = "C27001I_004",
                                          no_ins_19_64 = "C27001I_007",
                                          no_ins_65_up = "C27001I_010",
                                          households = "B11001I_001",
                                          family_hhs = "B11001I_002",
                                          married_couple_hhs = "B11001I_003",
                                          living_alone = "B11001I_008", 
                                          multigenerational = "B10051I_001",
                                          inc_und_10k = "B19001I_002",
                                          inc_15k = "B19001I_003",
                                          inc_20k = "B19001I_004",
                                          inc_25k = "B19001I_005",
                                          inc_30k = "B19001I_006",
                                          inc_35k = "B19001I_007",
                                          inc_40k = "B19001I_008",
                                          inc_45k = "B19001I_009",
                                          inc_50k = "B19001I_010",
                                          inc_60k = "B19001I_011",
                                          inc_75k = "B19001I_012",
                                          inc_100k = "B19001I_013",
                                          inc_125k = "B19001I_014",
                                          inc_150k = "B19001I_015",
                                          inc_200k = "B19001I_016",
                                          inc_over_200k = "B19001I_017",
                                          agg_income = "B19025I_001",
                                          owner_occupied = "B25003I_002",
                                          pplxroom_morethan1 = "B25014I_003",
                                          foodstamps = "B22005I_002",
                                          tot_educ_25up = "C15002I_001",
                                          lesshigh_m = "C15002I_003",
                                          lesshigh_f = "C15002I_008",
                                          highschool_m = "C15002I_004",
                                          highschool_f = "C15002I_009",
                                          somecollege_m = "C15002I_005",
                                          somecollege_f = "C15002I_010",
                                          bahigher_m = "C15002I_006",
                                          bahigher_f = "C15002I_011",
                                          occupation_tot_16_up = "C24010I_001",
                                          service_occ_male = "C24010I_004",
                                          production_occ_male = "C24010I_007",
                                          service_occ_female = "C24010I_010",
                                          production_occ_female = "C24010I_013",
                                          tot_pop_occupied_housing = "B25032I_001",
                                          units_oned = "B25032I_002",
                                          units_onea = "B25032I_003",
                                          units_two = "B25032I_004",
                                          units_three_four = "B25032I_005"))
        
        #### drop non NYC block groups
        acs_white <- mutate(acs_white, county = substr(GEOID,1,5)) %>%
                filter(county %in% NYC_counties) %>%
                select(c("GEOID", "variable", "estimate"))
        acs_black <- mutate(acs_black, county = substr(GEOID,1,5)) %>%
                filter(county %in% NYC_counties) %>%
                select(c("GEOID", "variable", "estimate"))
        acs_asian <- mutate(acs_asian, county = substr(GEOID,1,5)) %>%
                filter(county %in% NYC_counties) %>%
                select(c("GEOID", "variable", "estimate"))
        acs_hisp <- mutate(acs_hisp, county = substr(GEOID,1,5)) %>%
                filter(county %in% NYC_counties) %>%
                select(c("GEOID", "variable", "estimate"))
        
        #### spread and organize columns
        acs_white <- pivot_wider(acs_white, names_from = "variable", values_from = "estimate")
        acs_black <- pivot_wider(acs_black, names_from = "variable", values_from = "estimate")
        acs_asian <- pivot_wider(acs_asian, names_from = "variable", values_from = "estimate")
        acs_hisp <- pivot_wider(acs_hisp, names_from = "variable", values_from = "estimate")
        
        #### make percentage variables
        acs_white <- transmute(acs_white, 
                               tract = GEOID,
                               race = "white",
                               population = pop_white_nonhisp,
                               under17 = (male_und_5+male_9+male_14+male_17+fem_und_5+fem_9+fem_14+fem_17)/pop_white_nonhisp,
                               age18_44 = (male_19+male_24+male_29+male_34+male_44+fem_19+fem_24+fem_29+fem_34+fem_44)/pop_white_nonhisp,
                               age45_54 = (male_54+fem_54)/pop_white_nonhisp,
                               age55_64 = (male_64+fem_64)/pop_white_nonhisp,
                               age65_74 = (male_74+fem_74)/pop_white_nonhisp,
                               age75up = (male_84+male_85up+fem_84+fem_85up)/pop_white_nonhisp,
                               uninsured = (no_ins_18und+no_ins_19_64+no_ins_65_up)/pop_white_nonhisp,
                               households = households,
                               family_hhs = family_hhs/households,
                               married_couple_hhs = married_couple_hhs/households,
                               livingalone = living_alone/households,
                               multigenerational = multigenerational/households,
                               owner_occupied = owner_occupied/households,
                               more_people_than_rooms = pplxroom_morethan1/households,
                               foodstamps = foodstamps/households, 
                               ed_hs = (highschool_m + highschool_f)/tot_educ_25up,
                               ed_somecollege = (somecollege_m + somecollege_f)/tot_educ_25up,
                               ed_bach_up = (bahigher_m + bahigher_f)/tot_educ_25up,
                               ed_none = (lesshigh_m + lesshigh_f)/tot_educ_25up, 
                               service_occupation = (service_occ_male + service_occ_female)/occupation_tot_16_up,
                               production_occupation = (production_occ_male + production_occ_female)/occupation_tot_16_up,
                               inc_und_10k = inc_und_10k/households,
                               inc10_20k = (inc_15k + inc_20k)/households,
                               inc20_30k = (inc_25k + inc_30k)/households,
                               inc30_40k = (inc_35k + inc_40k)/households,
                               inc40_50k = (inc_45k + inc_50k)/households,
                               inc50_75k = (inc_60k + inc_75k)/households,
                               inc75_100k = (inc_100k)/households,
                               inc100_200k = (inc_125k + inc_150k + inc_200k)/households,
                               inc200_up = inc_over_200k/households,
                               inc_percap = agg_income/pop_white_nonhisp,
                               bldg_five_plus_units = (units_onea + units_oned + units_two + units_three_four)/tot_pop_occupied_housing,
                               institutionalized = (pop_white_nonhisp - tot_pop_occupied_housing)/pop_white_nonhisp)
        
        ####### black #########
        acs_black <- transmute(acs_black, 
                               tract = GEOID,
                               race = "black",
                               population = pop_black,
                               under17 = (male_und_5+male_9+male_14+male_17+fem_und_5+fem_9+fem_14+fem_17)/pop_black,
                               age18_44 = (male_19+male_24+male_29+male_34+male_44+fem_19+fem_24+fem_29+fem_34+fem_44)/pop_black,
                               age45_54 = (male_54+fem_54)/pop_black,
                               age55_64 = (male_64+fem_64)/pop_black,
                               age65_74 = (male_74+fem_74)/pop_black,
                               age75up = (male_84+male_85up+fem_84+fem_85up)/pop_black,
                               uninsured = (no_ins_18und+no_ins_19_64+no_ins_65_up)/pop_black,
                               households = households,
                               family_hhs = family_hhs/households,
                               married_couple_hhs = married_couple_hhs/households,
                               livingalone = living_alone/households,
                               multigenerational = multigenerational/households,
                               owner_occupied = owner_occupied/households,
                               more_people_than_rooms = pplxroom_morethan1/households,
                               foodstamps = foodstamps/households, 
                               ed_hs = (highschool_m + highschool_f)/tot_educ_25up,
                               ed_somecollege = (somecollege_m + somecollege_f)/tot_educ_25up,
                               ed_bach_up = (bahigher_m + bahigher_f)/tot_educ_25up,
                               ed_none = (lesshigh_m + lesshigh_f)/tot_educ_25up, 
                               service_occupation = (service_occ_male + service_occ_female)/occupation_tot_16_up,
                               production_occupation = (production_occ_male + production_occ_female)/occupation_tot_16_up,
                               inc_und_10k = inc_und_10k/households,
                               inc10_20k = (inc_15k + inc_20k)/households,
                               inc20_30k = (inc_25k + inc_30k)/households,
                               inc30_40k = (inc_35k + inc_40k)/households,
                               inc40_50k = (inc_45k + inc_50k)/households,
                               inc50_75k = (inc_60k + inc_75k)/households,
                               inc75_100k = (inc_100k)/households,
                               inc100_200k = (inc_125k + inc_150k + inc_200k)/households,
                               inc200_up = inc_over_200k/households,
                               inc_percap = agg_income/pop_black,
                               bldg_five_plus_units = (units_onea + units_oned + units_two + units_three_four)/tot_pop_occupied_housing,
                               institutionalized = (pop_black - tot_pop_occupied_housing)/pop_black)
        ####### asian #########
        acs_asian <- transmute(acs_asian, 
                               tract = GEOID,
                               race = "asian",
                               population = pop_asian,
                               under17 = (male_und_5+male_9+male_14+male_17+fem_und_5+fem_9+fem_14+fem_17)/pop_asian,
                               age18_44 = (male_19+male_24+male_29+male_34+male_44+fem_19+fem_24+fem_29+fem_34+fem_44)/pop_asian,
                               age45_54 = (male_54+fem_54)/pop_asian,
                               age55_64 = (male_64+fem_64)/pop_asian,
                               age65_74 = (male_74+fem_74)/pop_asian,
                               age75up = (male_84+male_85up+fem_84+fem_85up)/pop_asian,
                               uninsured = (no_ins_18und+no_ins_19_64+no_ins_65_up)/pop_asian,
                               households = households,
                               family_hhs = family_hhs/households,
                               married_couple_hhs = married_couple_hhs/households,
                               livingalone = living_alone/households,
                               multigenerational = multigenerational/households,
                               owner_occupied = owner_occupied/households,
                               more_people_than_rooms = pplxroom_morethan1/households,
                               foodstamps = foodstamps/households, 
                               ed_hs = (highschool_m + highschool_f)/tot_educ_25up,
                               ed_somecollege = (somecollege_m + somecollege_f)/tot_educ_25up,
                               ed_bach_up = (bahigher_m + bahigher_f)/tot_educ_25up,
                               ed_none = (lesshigh_m + lesshigh_f)/tot_educ_25up, 
                               service_occupation = (service_occ_male + service_occ_female)/occupation_tot_16_up,
                               production_occupation = (production_occ_male + production_occ_female)/occupation_tot_16_up,
                               inc_und_10k = inc_und_10k/households,
                               inc10_20k = (inc_15k + inc_20k)/households,
                               inc20_30k = (inc_25k + inc_30k)/households,
                               inc30_40k = (inc_35k + inc_40k)/households,
                               inc40_50k = (inc_45k + inc_50k)/households,
                               inc50_75k = (inc_60k + inc_75k)/households,
                               inc75_100k = (inc_100k)/households,
                               inc100_200k = (inc_125k + inc_150k + inc_200k)/households,
                               inc200_up = inc_over_200k/households,
                               inc_percap = agg_income/pop_asian,
                               bldg_five_plus_units = (units_onea + units_oned + units_two + units_three_four)/tot_pop_occupied_housing,
                               institutionalized = (pop_asian - tot_pop_occupied_housing)/pop_asian)
        ####### hispanic ###############
        acs_hisp <- transmute(acs_hisp, 
                              tract = GEOID,
                              race = "hispanic",
                              population = pop_hisp,
                              under17 = (male_und_5+male_9+male_14+male_17+fem_und_5+fem_9+fem_14+fem_17)/pop_hisp,
                              age18_44 = (male_19+male_24+male_29+male_34+male_44+fem_19+fem_24+fem_29+fem_34+fem_44)/pop_hisp,
                              age45_54 = (male_54+fem_54)/pop_hisp,
                              age55_64 = (male_64+fem_64)/pop_hisp,
                              age65_74 = (male_74+fem_74)/pop_hisp,
                              age75up = (male_84+male_85up+fem_84+fem_85up)/pop_hisp,
                              uninsured = (no_ins_18und+no_ins_19_64+no_ins_65_up)/pop_hisp,
                              households = households,
                              family_hhs = family_hhs/households,
                              married_couple_hhs = married_couple_hhs/households,
                              livingalone = living_alone/households,
                              multigenerational = multigenerational/households,
                              owner_occupied = owner_occupied/households,
                              more_people_than_rooms = pplxroom_morethan1/households,
                              foodstamps = foodstamps/households, 
                              ed_hs = (highschool_m + highschool_f)/tot_educ_25up,
                              ed_somecollege = (somecollege_m + somecollege_f)/tot_educ_25up,
                              ed_bach_up = (bahigher_m + bahigher_f)/tot_educ_25up,
                              ed_none = (lesshigh_m + lesshigh_f)/tot_educ_25up, 
                              service_occupation = (service_occ_male + service_occ_female)/occupation_tot_16_up,
                              production_occupation = (production_occ_male + production_occ_female)/occupation_tot_16_up,
                              inc_und_10k = inc_und_10k/households,
                              inc10_20k = (inc_15k + inc_20k)/households,
                              inc20_30k = (inc_25k + inc_30k)/households,
                              inc30_40k = (inc_35k + inc_40k)/households,
                              inc40_50k = (inc_45k + inc_50k)/households,
                              inc50_75k = (inc_60k + inc_75k)/households,
                              inc75_100k = (inc_100k)/households,
                              inc100_200k = (inc_125k + inc_150k + inc_200k)/households,
                              inc200_up = inc_over_200k/households,
                              inc_percap = agg_income/pop_hisp,
                              bldg_five_plus_units = (units_onea + units_oned + units_two + units_three_four)/tot_pop_occupied_housing,
                              institutionalized = (pop_hisp - tot_pop_occupied_housing)/pop_hisp)
        
        ### write files
        acs_full <- rbind(acs_white, acs_black, acs_asian, acs_hisp) %>%
                pivot_longer(cols = under17:age75up, names_to = "age", values_to = "pop_perc") %>%
                mutate(population = population*pop_perc,
                       strat_id = paste(tract, substr(race,1,1), substr(age, nchar(age) - 1, nchar(age)), sep = "_"))
        
        #### tract level stats, all races:
        
        acs <- get_acs(state = "NY", geography = "tract", year = i,
                       variables = c(population = "B01003_001",
                                     white_nonhisp = "B03002_003",
                                     black_nonhisp = "B03002_004",
                                     asian = "B03002_006",
                                     hispanic = "B03002_012",
                                     male_und_5 = "B01001_003",
                                     male_9 = "B01001_004",
                                     male_14 = "B01001_005",
                                     male_17 = "B01001_006",
                                     male_19 = "B01001_007",
                                     male_24 = "B01001_008",
                                     male_29 = "B01001_009",
                                     male_34 = "B01001_010",
                                     male_44 = "B01001_011",
                                     male_54 = "B01001_012",
                                     male_64 = "B01001_013",
                                     male_74 = "B01001_014",
                                     male_84 = "B01001_015",
                                     male_85up = "B01001_016",
                                     fem_und_5 = "B01001_018",
                                     fem_9 = "B01001_019",
                                     fem_14 = "B01001_020",
                                     fem_17 = "B01001_021",
                                     fem_19 = "B01001_022",
                                     fem_24 = "B01001_023",
                                     fem_29 = "B01001_024",
                                     fem_34 = "B01001_025",
                                     fem_44 = "B01001_026",
                                     fem_54 = "B01001_027",
                                     fem_64 = "B01001_028",
                                     fem_74 = "B01001_029",
                                     fem_84 = "B01001_030",
                                     fem_85up = "B01001_031",
                                     no_ins_und_18 = "B27010_017",
                                     no_ins_18_34 = "B27010_033",
                                     no_ins_35_64 = "B27010_050",
                                     no_ins_65_up = "B27010_066",
                                     households = "B11001_001",
                                     family_hhs = "B11001_002",
                                     married_couple_hhs = "B11001_003",
                                     one_person = "B11016_010",
                                     two_people_fam = "B11016_003",
                                     two_people_non_fam = "B11016_011",
                                     three_people_fam = "B11016_004",
                                     three_people_non_fam = "B11016_012",
                                     inc_und_10k = "B19001_002",
                                     inc_15k = "B19001_003",
                                     inc_20k = "B19001_004",
                                     inc_25k = "B19001_005",
                                     inc_30k = "B19001_006",
                                     inc_35k = "B19001_007",
                                     inc_40k = "B19001_008",
                                     inc_45k = "B19001_009",
                                     inc_50k = "B19001_010",
                                     inc_60k = "B19001_011",
                                     inc_75k = "B19001_012",
                                     inc_100k = "B19001_013",
                                     inc_125k = "B19001_014",
                                     inc_150k = "B19001_015",
                                     inc_200k = "B19001_016",
                                     inc_over_200k = "B19001_017",
                                     agg_income = "B19025_001",
                                     owner_occupied = "B25003_002",
                                     median_tenure = "B25039_001",
                                     tenure_1yr = "B25038_003",
                                     tenure_3yr = "B25038_004",
                                     tenure_8yr = "B25038_005",
                                     pplxroom_1_1.5o = "B25014_005",
                                     pplxroom_1.5_2o = "B25014_006",
                                     pplxroom_2pluso = "B25014_007",
                                     pplxroom_1_1.5r = "B25014_011",
                                     pplxroom_1.5_2r = "B25014_012",
                                     pplxroom_2plusr = "B25014_013",
                                     public_asst = "B19057_002",
                                     tot_transp = "B08301_001",
                                     public_transp = "B08301_010",
                                     tot_educ_25up = "B15003_001",
                                     highschool = "B15003_017",
                                     ged = "B15003_018",
                                     some_college_less1yr = "B15003_019",
                                     some_college = "B15003_020",
                                     associate = "B15003_021",
                                     bachelor = "B15003_022",
                                     master = "B15003_023",
                                     professional = "B15003_024",
                                     doctorate = "B15003_025",
                                     tot_rent2inc = "B25070_001",
                                     no_rent2inc = "B25070_011",
                                     rent2inc_50_up = "B25070_010",
                                     rent2inc_40_50 = "B25070_009",
                                     rent2inc_35_40 = "B25070_008",
                                     rent2inc_30_35 = "B25070_007",
                                     occupation_tot_16_up = "C24010_001",
                                     service_occ_male = "C24010_019",
                                     production_occ_male = "C24010_034",
                                     service_occ_female = "C24010_055",
                                     production_occ_female = "C24010_070",
                                     tot_pop_occupied_housing = "B25033_001",
                                     bldg_units_5_up_own = "B25033_005",
                                     bldg_units_5_up_rent = "B25033_011",
                                     gini = "B19083_001"))
        #### drop non NYC block groups
        acs <- acs %>% mutate(county = substr(GEOID,1,5)) %>%
                filter(county %in% NYC_counties) %>%
                select(c("GEOID", "variable", "estimate"))
        #### spread and organize columns
        acs <- pivot_wider(acs, names_from = "variable", values_from = "estimate")
        #### make percentage variables
        acs <- transmute(acs,
                         tract = GEOID,
                         tot_population = population,
                         white = white_nonhisp/population,
                         black = black_nonhisp/population,
                         asian = asian/population,
                         hispanic = hispanic/population,
                         other_nonhisp = 1 - (white + black + asian + hispanic),
                         tot_under17 = (male_und_5+male_9+male_14+male_17+fem_und_5+fem_9+fem_14+fem_17)/population,
                         tot_age18_44 = (male_19+male_24+male_29+male_34+male_44+fem_19+fem_24+fem_29+fem_34+fem_44)/population,
                         tot_age45_54 = (male_54+fem_54)/population,
                         tot_age55_64 = (male_64+fem_64)/population,
                         tot_age65_74 = (male_74+fem_74)/population,
                         tot_age75up = (male_84+male_85up+fem_84+fem_85up)/population,
                         median_tenure = 2018 - median_tenure,
                         tenure_8yrs_under = (tenure_1yr + tenure_3yr + tenure_8yr)/population,
                         tot_uninsured = (no_ins_und_18+no_ins_18_34+no_ins_35_64+no_ins_65_up)/population,
                         tot_households = households,
                         tot_family_hhs = family_hhs/households,
                         tot_married_couple_hhs = married_couple_hhs/households,
                         tot_big_hhs_4_plus = 1 - ((one_person + two_people_fam + two_people_non_fam + three_people_fam + three_people_non_fam)/households),
                         tot_owner_occupied = owner_occupied/households,
                         tot_pplxroom_1_2 = (pplxroom_1_1.5r + pplxroom_1_1.5o + pplxroom_1.5_2r + pplxroom_1.5_2o)/households,
                         tot_pplxroom_2plus = (pplxroom_2plusr + pplxroom_2pluso)/households,
                         tot_pplxroom_less1 = 1 - tot_pplxroom_1_2 - tot_pplxroom_2plus,
                         tot_public_assistance = public_asst/households,
                         tot_public_transport = public_transp/tot_transp,
                         tot_ed_hs = (highschool + ged)/tot_educ_25up,
                         tot_ed_assc_somecollege = (some_college_less1yr + some_college + associate)/tot_educ_25up,
                         tot_ed_bach = bachelor/tot_educ_25up,
                         tot_ed_higher = (master + professional + doctorate)/tot_educ_25up,
                         tot_ed_none = 1 - (tot_ed_hs + tot_ed_assc_somecollege + tot_ed_bach + tot_ed_higher),
                         tot_rent2income50percent = rent2inc_50_up/(tot_rent2inc - no_rent2inc),
                         tot_rent2income30_50percent = (rent2inc_30_35 + rent2inc_35_40 + rent2inc_40_50)/(tot_rent2inc - no_rent2inc),
                         tot_service_occupation = (service_occ_male + service_occ_female)/occupation_tot_16_up,
                         tot_production_occupation = (production_occ_male + production_occ_female)/occupation_tot_16_up,
                         tot_inc_und_10k = inc_und_10k/households,
                         tot_inc10_20k = (inc_15k + inc_20k)/households,
                         tot_inc20_30k = (inc_25k + inc_30k)/households,
                         tot_inc30_40k = (inc_35k + inc_40k)/households,
                         tot_inc40_50k = (inc_45k + inc_50k)/households,
                         tot_inc50_75k = (inc_60k + inc_75k)/households,
                         tot_inc75_100k = (inc_100k)/households,
                         tot_inc100_200k = (inc_125k + inc_150k + inc_200k)/households,
                         tot_inc200_up = inc_over_200k/households,
                         tot_inc_percap = agg_income/population,
                         tot_bldg_five_plus_units = (bldg_units_5_up_own + bldg_units_5_up_rent)/tot_pop_occupied_housing,
                         tot_institutionalized = (population - tot_pop_occupied_housing)/population,
                         gini = gini)
        
        name <- paste("acs", i, "_nyc_tract_clean_byrace.csv", sep = "")
        write.csv(acs_full, name)
        
        name2 <- paste("acs", i, "_nyc_tract_clean.csv", sep = "")
        write.csv(acs, name2)
}
        

